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SUMMARY 

/I theoretical flutter analysis has been made that 
utilizes the same simplified model employed in the 
vibrational investigation presented in TN D-087 in 
which chordwise and spanwise flexibilities of the 
actual wing are duplicated . This idealized model 
subjected to a supersonic flow is investigated to 
predict the flutter behavior and establish a flutter 
boundary . Additional effects such as aerodynamic 
damping and variations in cross-sectional curvature 
are determined and their effect on the predicted 
flutter boundary ascertained . A comparison of the 
predicted flutter boundary with results obtained 
experimentally for multiweb wing structures shows 
very good agreement. 

INTRODUCTION 

In an experimental investigation conducted by 
the Langley Structures Research Division into the 
dynamic flutter failures of multiweb wing struc- 
tures, chordwise or camber flut ter was found to be 
a primary mode of failure when insufficient bulk- 
heads were used. This flutter consisted of small- 
amplitude oscillations with five or less spanwise 
node lines across the chord. Tn the investigation a 
scries of multiwch wing structures were tested in 
the preflight jet of the NASA Wallops Station at a 
Mach number of 2 as discussed in references 1, 2 > 
and 3. Additional tests have also been made in 
the Langley 9- by 6-foot thermal structures tunnel 
at a Mach number of 3, and these flutter data are 
included also. These multiweb models had 
5-percent -thick circular-arc airfoil sections con- 
taining six webs and a solid leading and trailing 
edge. Tn attempts to predict this flutter behavior 
theoretically, modal approaches using the lowest 
five experimental mode shapes have not been very 
successful. Vibration behavior of such structures, 


however, has been successfully predicted theo- 
retically (for example, ref. 4) by the use of an 
analysis of a simplified model which took cogni- 
zance of the cross-sectional flexibilities. Hence, 
the same idealized model subjected to a supersonic 
flow lias been investigated in an attempt to obtain 
a method of predicting the flutter behavior of the 
multiweb wings. D. J. Johns, in a note on the 
influence of panel deformations on wing flutter 
(ref. 5), also recognized the significance of cou- 
pling between cross-sect ion al and spanwise defor- 
mation modes. This analysis, however, was 
applicable only to wings with rectangular cross 
sections and consequently did not include any 
effects of curvature of the cross section. 

The purpose of the present paper is to present 
the results of a flutter analysis of the idealized 
structure and to show how these results can be 
used to predict the flutter behavior of multiweb 
and monocoque low-aspect-ratio wings of circular- 
arc cross section. The idealized structure consists 
of a typical cross section of the wing mounted 
on springs and lienee embodies the principal 
mechanisms of the cross-sectional and spanwise 
flexibilities of the actual wing. A modified linear 
piston theory is used to represent the aerodynamic 
forces created by supersonic flow. The analysis 
consists in solving the basic equations by means of 
a Galerkin procedure using a modal approximation 
of three and four terms but neglecting the effect of 
aerodynamic damping. The accuracy of this 
Galerkin procedure is determined in two ways — 
first, by the apparent convergence of the results 
of the three- and four-term solutions and second, 
by a comparison of the flutter boundaries for the 
case of the idealized structure with infinitely stiff 
spring supports (pin-ended) at the leading and 
trailing edges as calculated by the Galerkin method 
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with those obtained from an exact analysis. In 
addition, the effects of damping are evaluated by a 
two-mode Galerkin solution that includes the 
effects of aerodynamic damping. Finally, the 
results of the analysis are used to predict suc- 
cessfully the flutter boundary of the multiweb 
wings used in the experimental investigations. 
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SYMBOLS 
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E 
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h 
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j, 1 
k 

M 
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<L 

t 

U m 

iv 

w 

X, Ih z 

a 

7, 5, 9, 0 


A 

A„ 

X 

X 

Hr 

P 

Pco 

T 

*=£-1 1 


chord length of beam 
modulus of elasticity 
distance from /-axis to median line of 
skin 

value of h at /~0 

area moment of inertia about median 
line of skin, per unit length, ff ] 2 
area polar moment of inertia per unit 


length, 


E 


integers 
spring constant 
Mach number 

roots of auxiliary equation (see ap- 


pendix) 


1 


dynamic pressure, - p^u 

skin thickness 
free-stream velocity 
amplitude of deflection in ^-direction 
(positive upward) 
deflection in ^-direction, Re(w fwT ) 
Cartesian coordinates (fig. 1) 

frequency coefficient, 

\I6 hi J 

assumed coefficients of the roots of 
the exact equation (see eqs. (A9) of 
appendix) 

dynamic-pressure parameter, (Q 
critical value of A 

amplitude of normal midplane force 
in the beam (positive in tension) 
normal midplane force, Re(Xc ?wr ) 
Poisson's ratio 


mass density of cover 

mass density of air (free stream) 

time 


« natural frequency 

Subscripts: 


A ant isymmetrical 

i f m, n any integer 

S symmetrical 

T torsional 

x, y, 2 , t variable with respect to which primary 

symbol is differentiated 

Primes and superscript IV on symbols denote dif- 
ferentiation with respect to £. 


ANALYSIS 

IDEALIZED STRUCTURE 

A rigorous approach to the flutter analysis of t lie 
multi web wing shown in figure 1(a) would, of 
necessity, involve some type of shell analysis and 
sophisticated aerodynamics. A more convenient 
though less rigorous approach is to investigate an 
idealized structure identical to that used in the 
vibration analysis of reference 4. The only struc- 
tures considered are those that have relatively uni- 
form properties in the span vise direction with a 
doubly symmetrical curved air foil cross section as 
shown in figure 1(a). The wing or beam is as- 
sumed to be of either* pure monocoque or multi web 
cross section with no or relatively few internal ribs. 

The idealized structure consists of a typical cross 
section of the wing mounted on identical elastic 
supports as shown in figure 1(b). The cross sec- 
tion consists of two curved uniform beams at- 
tached at the ends in such a manner as to maintain 
the angle between the beams at each end. The 
rigid-body motion of the idealized structure de- 
pends on the stiffness of the elastic* supports which 
simulate the beam bending and torsional stiffness 
of the wing. The depth of the individual beams 
in the idealized structure is taken to be the same 
as the thickness t of the covers in the prototype. 
The curvature of the covers of the wings consid- 
ered is such as can be represented by the para- 
bolic* equation k=ho (] — £ 2 ) where ho is the maxi- 
mum value of h and £ is the non dimensional co- 
ordinate 2z/c. 

The webs of the multi web structure are assumed 
to contribute nothing to the camber stiffness. 
They do, however, tend to prevent any relative 
motion between the top and bottom covers. Thus, 
the additional assumption is made (fiat the deflec- 
tions of the top and bottom beams in the idealized 
structure are the same. Each beam of the ideal- 
ized structure is permitted to carry midplane forces 
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(a) Prototype, 
(b) Idealization. 


Fku’re I.- Low-aspect-ratio multiweb wing. 

in addition to transverse shear forces and bending 
moments and is considered to be acting in a man- 
ner consistent with elementary theory (i.e., plane 
sections remain plane). The air loads to which 
the idealized struct urc are subjected are taken as 
those given by modified linear piston theory. 


where 

U m free-stream velocity 

q dynamic pressure, 

p compressibility factor, yM 2 — 1 

and the subscripts x and r denote differentiation 
with respect to x and r. The midplane force X, 
as in reference 4, was assumed constant over the 
length of the beam (for any finite r) and an inte- 
gral expression for X was arrived at by the as- 
sumption that no overall shortening or lengthening 
of the beams was possible. Thus, 

- Ft r c/2 

htWjdx (5) 

C — c/2 

Note that when w is an tisym metrical, the inte- 
grand of equation (5) is an an Asymmetrical 
function about the midchord (3=0); consequently, 
its integral from —c/2 to c/2 is zero and no mid- 
plane forces are present. This is not true for the 
symmetrical case. 

Tire differential equation and associated boun- 
dary conditions are similar to equations (9) to (13) 
of reference 4 except that terms are included on 
the right-hand side of equation (1) to represent 
the air loads given by modified linear piston 
theory (modified in the sense that 1/M is replaced 
by 1/0). The first term on the right-hand side 
of the mod ified-pist on- theory equation (1) is a 
static-loading term and the second an aero- 
dynamic damping term. The solution of equa- 
tion (1) can be written in the form 

ic = Tte f?r(/)c lwr | (0) 
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-0 






£--i 


=0 


(9) 

( 10 ) 

(11) 


where the primes and the superscript IV denote 

2 x 

differentiation with respect to £=— ? and A is a 

c 

dynamic -pressure parameter defined as 

-IKO' <> 2 > 

A frequency parameter 7? 4 is defined as 


n 4 3 qc 4 . . 

B ~2E0i 3 U m lu a 

4 3c 4 pw a 


4Et 2 


(1.3) 


(H) 


equation and boundary conditions are identical to 
that given in reference 7 for an infinite-aspcct- 

mtio plate (with the exception of a factor of — * 2 

1 fX 

resulting from the use of the plate stiffness D in 
the panel analysis and El in the beam analysis). 
If the aerodynamic damping is neglected, the 
last term in equation (1) is omitted and the ex- 
pression for B 4 given by equation (13) reduces to 
— a 4 . The solution to this differential equation 
and boundary conditions is discussed in the ap- 
pendix. In the appendix it is shown that the de- 
flection w can be written as 

™sh H+C, sinh <t>\p) 

+e- e +(C 3 cos ty+C A sin Sf) (16) 

where the relations between <f> , 8 } a 4 , and^Tare given 
by the following equations: 


and 



(17) 

2Et f ' 4Eth n ( 

X=-^- j ^ h w f/£=- c2 



(IS) 

EXACT SOLUTION 

(15) 

“’-it?-*" 

(19) 

An exact solution to equation (1) in which aero- In addition 

dynamic damping has been neglected is obtained 

A 1/3 ^ 


in this section. The approach is 

similar to that 

2“< 0 <° l 

(20) 

used in reference 6, which presents an exact solu- 
tion for panel flutter of a simply supported plate. 

H>l*l J 



Although solutions of the flutter of the idealized 
structure considered in this section are possible for 
any value of spring constant k, only the special 
case of infinite k is presented. This approach re- 
sults in considerable simplification in boundary 
conditions and consequent numerical procedures; 
for example, the solutions are not dependent on 
values of cjt or E. Although considerably simpli- 
fied, the solution still provides a means of evalu- 
ating the effect of curvature and its associated 
normal midplane force X on flutter behavior. 
These are the principal effects that distinguish the 
present analysis from typical wing-flutter or panel- 
flu tier studies. If no initial curvature is present 
as in the case of a straight beam, h 0 is considered 
zero and consequently X is zero, and the differential 


Substitution of equation (1G) into the boundary 
conditions (eqs. (A4) and (A />)) results in the 
characteristic equation presented in table I. 
Although as written in table I the determinant is 

a function of 0, <j> } 5, a, A, and ^ > the parameters 

4>, 5, and a can in turn be written as functions or 
0 and A, (Sec eqs. (17) to (19).) Consequently 
the determinant can be considered as a function of 


value of A and 


only 0, A, and * For a given 

( j^ fhe appropriate value of 6 can be determined 

from the characteristic equation. The corre- 
sponding value of a can then be computed from 
equation (19). 



CAMBER FLUTTER OF MONOCOQUE AND MULTIWEB WINGS 


5 


TABLE I— FLUTTER DETERMINANT OF EXACT SOLUTION, WHEN AND AERODYNAMIC 

DAMPING IS NEGLECTED 


whore 


1 +alh 
(0 2 + * 2 ) 

\<> 29 cosh 20—1 

fh 


aJh 

200 

e 26 si nh 2 0 

Th 


1 


IU = 


u,=- 


ih=- 


c n ^eosh 2 0 — ^ si nh 2<^>^ — ■ 1 

Kf 

^ si nil 2<^> — ^ cosh 20^+^ 

KF 

2 * ^ — cos 26-i-| sin 26^ + 1 


/f 4 = - 


■*(■ 


1 + 


-sin 25 — ^ cos 25 ) + ^ 


1 + 


(i ) 1 


In figure 2 a plot of the dynamic-pressure param- 
eter A as a function of the frequency parameter a 


is shown by the solid curves for 


(0 


= 0. For 


discrete values of A an infinite set of values of a 
is shown to exist. When A equals zero, the values 
of a correspond to the natural vibration frequen- 
cies determined in reference 4 for a monocoque 
beam in a vacuum. As A increases from zero, 
these natural frequencies change until a value of 
A is reached in which two of the roots of a coalesce. 
For values of A greater than this value two of the 
values of a become complex and hence, as can be 
seen from equation (14), to will be complex. Since 
o) must remain real to give stable oscillatory mo- 
tion, this value of A is the critical value A cr and 
thus defines the flutter boundary. In figure 3, 

A c , is plotted as a function of ( from 0 to 0.25. 
//V V * o/ 

Note that ^ J for all practical configurations is 

much less than unity. For example, in an extreme 
case where the skin thickness is equal to half the 
value of h Q , this ratio is only 0.25. In the more 
practical type of construction of multiweb wings 
such as those presented in references 1, 2, and 3 


1 + aII, 
(0 2 ~5 2 ) 

e~ 28 cos 26—1 
Th 


aTh 
-205 

c~ 2e sin 25] 

Th 


J/. = C 2*[(02 j 02 ) (, 0 ^ 20 + 200 sinh 2 0 ] 


// 6 =c 2 *[( 0 2 + 0 2 ) sinh 20 -f 200 cosh 20 ] 


2 , [(0 2 -5 2 ) cos 25 + 205 sin 26] 


Th 2 ^( 02 - 52 ) s i n 25-206 cos 26] 


A 2 

: ic>0 2 ‘ 


• is 4 



Fiocre 2. — Effect of dynamic-pressure parameter on pin- 
ended (k= oo ) wing frequencies for = 0 * 

/3c 4 pto 2 V/L 3f/ / cV 

a ^\TEt 2 ) ; ‘ t) ' 

the ratio (^j ^ is in the range of 0.020. Con- 
sequently, within the range of practical design of 
multiwcl) and monocoque wings the ratio 
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(£) 2 


Figure 3. KfTect of variation in ^ on for exact 
solut ion. k= co . 


can bo assumed to approach zero without appre- 
ciably a flee ting the value of A rr . For tliis reason 
in the analysis performed by the Galerkin pro- 
cedure in the next section, a constant value of 



0 is used as tliis assumption results in con- 


siderable simplification in the numerical computa- 
tions. Although, as was previously discussed, 

^ ^ for all practical cases is much less than one, 


a unique situation occurs in the limiting case 
where the two curved beams merge into a straight 
beam. For this limiting case h 0 is assumed to be 


zero and 


(0 


tt can be shown that for this 


case the differential equation and boundary con- 
ditions (eqs. (I), (A4), and (Ao)) limit to the 
equations presented in reference 7 for the iufmitc- 
aspcct-rntio fiat panel. The value of A cr for 



is 37.90 and is indicated in figure 3 in 


order to compare the effects of cross-sectional 
curvature and its associated midplane forces with 
flat-panel results containing no in-plane forces. 
Comparing this value of A cr with that calculated 

for ^ ^ =0 (A cr ^-70.9) demonstrates that the 


existence of curvature can have appreciable effect 
on the flutter criterion. However, as can be seen 
from figure 3, once the curvature is present, varia- 
tions in the degree of curvature can have little 
effect on the final flutter boundary. 


SOLUTIONS BY THE GALERKIN PROCEDURE 

The previous solution, although exact, was for 
the special case of pin-ended constraint conditions 
where k= c° . Tu order to study the effect of 
span wise flexibility on the overall motion of the 
idealized structure finite values of the elastic 
supports must be used. Because of the resulting 
complexity in t lie final determinant which occurs 
in the case of an exact solution, a Galerkin pro- 
cedure using a theoretical modal approach is used 
to determine the flutter characteristics. As was 
done in the exact solution, the effects of aero- 
dynamic damping are neglected. The importance 
of this effect is discussed in a subsequent section. 

In order to perform a Galerkin solution by a 
modal approach, the choice of modes to be applied 
must satisfy both the geometric boundary con- 
ditions and the natural boundary conditions. 
Since the inclusion of the dynamic pressure q in 
equation (1) does not alter the boundary condi- 
tions associated with the elastic supports, the 
mode shapes and frequencies of the vibrational 
analysis of such structures derived in reference 4 
are used in the Galerkin procedure. That is, the 
deflection is represented by 

/ i 

w — (i-qWQ' s + 2 a»w„ t s + Xj b m w mtA ( 21 ) 

n = 1 m =(] 


where Uh.s, and w, Uf A are the symmetric 

and antisymmetric vibration modes. Note that 
the coefficients a 0 , a r(} and b m can be complex but 
WO'S, w nt s, and u\ n A are real. The mode shapes 
w r)% fi and w mtA and their associated frequency 
equations have been derived in reference 4 for all 
values of t/h ti . Examination of the equations 
given in reference 4 shows that a significant sim- 
plification results when f/A 0 approaches zero. 
Inasmuch as the exact solution has shown that 
this simplifying assumption would have no appre- 
ciable effect on the resulting flutter boundary for 
k=°°, tjh o is assumed to be zero, whereas eft 
remains finite. The frequency equation for the 
symmetrical case as obtained from reference 4 for 

O’ ■ 


cronies 


2 a«, 5 =tan tanh 


( 22 ) 


and the associated mode shape is written as 
... _ cos , cosh or*. <*£ 0 

w?n,S — __ I 1 " 

cos a ttiS COSll a„, s 


(2.3) 
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where 



For the anlisymmelrieal modes, as was pre- 
viously discussed, the midplane force X is zero, 
and the differential equation of vibration simply 
reduces to the elementary beam equation. Thus, 
the characteristic or frequency equation can be 
written as 


coth o',,, a co t .1 _ ‘3k / cY 

2 CCm, A 4 Kal , , A V / 

and the mode shape is 

^ sm (Xm., i smh 

sin a m , A ' sinh a m , A 

where 

/3 C + pa>«, a\ 1/4 
a,n - A \ 4 Et 2 ' ) 


(25) 


(2G) 


(27) 


It can be shown by the use of the general mode 


shapes given in reference 4 that regardless of the 
value of (j[ y, 1,10 natural vibration modes are 
orthogonal for all values of k. Consequently, 



w„,sWi,sd £= 0 

iv„,sWm, A dZ = 0 
w m , A w UA <lt = 0 


(i 9^rt) 


/(»= 0 , 1 , 2 , 

\ ((= 0 , 1 , 2 , 


• • 3) 

■ J) 
(28) 


/ (n=0, 1,2,... j) 
\(m=0, 1,2, ... 1) 


( 20 ) 


( i m ) 


f H= 0, 1,2,.../) 

\(m=0, 1,2,.../) 


(30) 


Note from equation (23) that for (yt) =0, the 

tip deflections of the cross section in the sym- 
metrical modes are zero and 



(n—l, 2, . . . j ) 


(31) 


and 


I W ”‘,^ W n,^h j V) m , A W n ,S<k ^( ?n _0, 1, 2, . .' /) ^ 

In addition, with the use of equations (23) and (2G) the following relationships are obtained 

j* wSLsd^tan 2 a H ,s— tiinh 2 (w=0, 1, 2, ... j) (33) 

J* W?«.a^=2+C0t 2 a m , a— colli 2 ot,n t (m=0, 1,2,.../) (34) 




, , OC„ S , , , (•/' / C\ 3 ( 

wm, au\ y/f = 4 — (««. ^ cottl «*. -i— b +Tc l 7 1 l ; v 4 m 

<X m , .4. — °^n, S ^ \t/ \«m,A i 


8«i.i 


' 2 AO&,S 


<x,i, s tan g w , A / (m = 0, 1 ,2, . . . 

ai.s / \ (^=0, 2, . . . 


0 

i) 


(35) 


The first mode, w 0tS , associated with a 0 in equa- 
tion (21) is the rigid-body translational mode 
while w {)tA consists of predominantly rigid-body 
torsional or pitching motion; that is, motion of 
the cross section as a rigid body mounted on 
spring supports. The- rigid-body torsional mode 
is always the lowest of the ant asymmetrical modes, 
while the rigid-body translational mode is not 
necessarily the lowest of the symmetric modes. 

6735S0 — 63 2 


(Sec ref. 4.) The rest of the modes associated 
with a n and b m are predominantly cross sectional 
in nature containing two or more nodal points. 

By application of the Galerkin procedures and 
use of the simplifying orthogonality relations given 
by equations (28) to (30) and equations (31) and 

(32) for ^“^=0, a set of equations for the 

coefficients a n and b m are found. The set of 
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equations Tor the coefficients a n and b m are: 


(«#. s— « 



'l£—0 


(a«, s 



s<ta- A 



(ai. .4 —a')b m \ id. + A fl„ f IP*, ,s'/| = 0 

*/ —I w=l J-l 


f («= 1, 2, . . . j) [ 
\(m-0, 1 , 2 , . . . I) f 

r(w=0J,2, .. . 1 ) 

\ (« = 1,2 


( 30 ) 


It can be seen from equations (3G) tlmt (at least in the absence of aerodynamic damping) the rigid-body 
symmetrical (translational) mode n\ s becomes uncoupled in the solution. Furthermore, it can be 
shown from equations (3G) that if the expansion for w in equation (21) contains either all symmetric 
or all antisymmetric modes, the resulting solutions do not yield an instability; lienee, flutter occurs 
only through the coupling of symmetric and antisymmetric modes. For a nontrivial solution, the deter- 
minant of the coefficients must be equal to zero. 

Three-term solution. — For a three-term approximation in which the rigid -body torsional mode and 
the first two cross-sectional modes are considered, the determinant of equations (3G) reduces to 


o 1 ) | «>!,*•'/{ 

— A j" Wi tA w' us 'IS 

— A ^ n\ A h 

Pi.Wf 


A j ( ii’i. aw[. 

(of*. -4 — Of 4 ) j* ^ w\ tA d% 

0 


= 0 

A ( ^ w 2 ' A w[.s<1i 

0 

(ai A —a 4 ) ^ ' 

Wlji; 

i 


which upon expansion gives 

(aj, *— a 4 )(af. a — a 4 ) (<4 a —a 4 ) |* > IP?, ^ w\. A d$ j* ipjj, A <% 

(a* aj, ^t) ( ^ W it A w[, ,sd) j* ^ IP?, a<H+ (a 4 — a4..0 ( |* | U’i, .1 IP,', .s^) j" ^ iPjL^/f 


(37) 


(38) 


The variation of the frequency parameter a 
for finite values of the dynamic-pressure parameter 
A can be obtained from equation (38) by choosing 
values of a and solving for A. Equation (38) also 
is a function of the wing material 77, chord- 
thickness ratio eft, and the spring constant k as 
can he seen from the values given for the integrals 
(eqs. (33) to (35)). Tn figure 2 a plot or A as a 
function of a is shown for the three-term Galcrkin 
approximation (by the long-short dashed curves) 


for 


■ ( 0 - 


0 and k 


As in the ease of the 


exact solution the critical value of A is the lowest 
value of A at which two of the frequency param- 
eters a coalesce. 

The critical value of A, expressed in the form of 


the often used panel 


flutter parameter 


t 

c\<i/ 



1/3 


is shown plotted in figure 4 as a function 


of k, for an aluminum-alloy (77— 10.GX 10* psi) 
wing of chord length c of 20 inches and a ratio eft 
of 312.5. This ratio eft was chosen as it corre- 
sponds to the ratio for the multiweb wings used 
in the experimental investigations. The effect of 
varying eft on the final flutter boundary will be 
discussed in a forthcoming section. The region 
above the computed curve represents the region 
of stability while the region below represents the 
region where flutter w ill occur. 

In order to convert t ho abscissa in figure 4 to a 
nondimensional form, the results of figure 4 were 
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Figure 4.-- Variation of the flutter parameter, --- with k for an aluminum-alloy wing (77 — 10.6X10° psi) with 

' (?)“-(£)“’ (O'-- 


4=312.5. J 
t c 


replotted as a function of the square of the 
frequency ratio cc T /co ltS in figure 5. The fre- 
quency co 7 > is t lie rigid-body torsional or pitching 
frequency of the idealized structure, whereas 
co l>s is the frequency of the first symmetrical 
camber mode. The uncoupling of the rigid-body 
transverse mode a motivated the choice of co T 
for use in the nondimensional frequency ratio. 
Thus, 

T\. f 


where 


I,P _9 k /rV 
1 < ) Kla\ s 4 $ \t/ 

(39) 

Cpt 


KrM O’ 

(40) 

r =Tl 

(41) 


and the approximate polar moment of inertia I p 


is given by 


I p ~I yi 


2 te 
= 12 


(42) 


Note from t lie curves presented in figure 5 that 
the torsional beam stiffness of the prototype as 
represented by co T is relatively unimportant for 
large values of co r and furthermore that flutter will 

occur even when <o T is infinite ^three-term solution, 
- ^ —0.348/)^. Tlic reversal of the curve at 

_2.i (or k— 145) for the three-term Galer- 

\u h s/ 


kin approximation is similar to that occurring in 
other solutions that use three-mode approxima- 
tions such as flat-panel modal solutions as shown 

in reference 7. For values of ( /aJ? N )>2.1 the 

\wi.s/ 
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critical flutter boundary is defined by the coales- 
cence of the rigid-body torsional (first an Asym- 
metrical) and the first symmetrical chordwisc 

mode while for } <"2.1 the first symmetric 

\wi ,s/ 

and second an Asymmetrical chordwisc modes 
coalesce. (See fig. 6(a).) As ( 0)T ) is further 

\" 1 , 5 / 

decreased, a point is reached at which the second 
antisymmetric chordwisc frequency becomes equal 
to the first symmetric chordwise frequency for 
zero airflow; at this point it appears as though 
flutter will occur at zero dynamic pressure. This 
phenomenon is quite common in flutter investiga- 
tions and is probably due to the neglect of aero- 
dynamic damping. 

Four-term solution. — In order to investigate 
the accuracy in the region of the curve reversal 
and also the convergence of the three-term approx- 
imation, a four-term Galerkin solution was per- 
formed. For a four-term approximation in which 


the fust two an Asymmetrical chordwise modes 
(one of which is predominantly rigid-body torsion) 
and the first two symmetrical chordwise modes 
(the rigid-body transverse mode having become 
uncoupled in the solution) are considered, the 
determinant of equations (36) upon expansion 
reduces to a quadratic in A 2 . Solving for the 
critical value of A by varying a results in two 
peaks, the minimum peak being A cr . 

The correlation between these approximate 
solutions and the exact solution for the specific 
case of pin-ended constraint conditions and 

GO * s s ^ own ’ n fi£ urc 2 * The agreement 

between the three- and four-term approximations 
and the exact solution— is- excellent, the critical 
flutter parameter being calculated as 70.5 for 
the three-term solution, 70.72 for the four-term 
solution, and 70.9 for the exact solution. In fad, 
the values of the curves as calculated for the first 
coalescing modes are so close that the curves 
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(L) 2 = o. £=312.5; F. IO.GX KF psi. 
\/?y/ t 


drawn through them cannot be easily distinguished 
one from another. The difference between the 
exact peaks and the approximate peaks becomes 
evident only at the higher coalescing modes. This 
agreement demonstrates the. degree of excellence 
in comparison between the approximate Galerkin 
procedures and the exact solution for the special 
case of infinite k. 

For the case of finite k the rapid convergence of 
the Galerkin procedures using only three and four 
modes, as is shown in figures 4 and 5, indicated 
that the inclusion of more terms in the approxi- 
mate solution was not necessary in order to obtain 


accurate results. 

Now that the accuracy of the Galerkin pro- 
cedure has been evaluated it is of interest to 
examine more closely the behavior of the partici- 
pating modes in the vicinity of the curve reversals 
shown in figures 4 and 5. In order to illustrate 
the change in coalescing modes, A has been plotted 
as a function of a for different values of k in 
figure 6(a) for the three-term approximation and 
in figure 6(b) for the four-term approximation. 
By following the sequence of plots in figure 6(a) 
for different values of k it can be seen that as the 
value of k varies the magnitude and shape of the 
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(fo) Four-term Oalerkin solution. 
Ftgt tie 6 . —Concluded, 


peak changes until different modes coalesce. As 
was previously discussed the reversal of the curve 
for the three-term approximation was to he 
expected because of t he cubic nature of the ex- 
panded determinant, as was demonstrated in flat- 
panel modal solutions. (See ref, 7.) The abrupt 
change in the curve for the four-term approxima- 
tion shown in figure 5 occurring at f - w -- j =1.93 

(k— 133) is somewhat surprising, however. This 
abrupt change for the four-term approximation is 
due to the type of modal coalescence illustrated 
in the sequence of plots in figure 0(b). Tn figure 
0(b), as in figure 0(a), the variance in k causes a 
change in the magnitude and shape of the two 
peaks until different modes coalesce. The coales- 


cence at values of ( -)]>]. 93 occurs between 

V^i ,s/ 

the rigid-body torsional (first an Asymmetrical) 
and the first symmetrical chordwisc mode, and 
between second symmetrical and second anti- 
symmetrical chordwisc mode. When the value 

of k decreases until -^<A.93, then the rigid-body 

w i,s 

torsional mode coalesces with the second sym- 
metrical chordwisc mode and beneath this curve 
the second an Asymmetrical and first symmetrical 
chordwisc modes coalesce, the smaller of these two 
peaks yielding A cr as always. 

Effects of variation in chord-thickness ratio. — 
Tn the previous Galerkin procedures the ratio of 
the chord to the skin thickness cjt was held 
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constant at 312.5, the value computed from the 
physical characteristics of the multiweb wings 
investigated in references 1 and 3. 

in order to show the effect of variance in the 


ratio eft on the flutter criterion, values of 


t/Epy* 

C \(£ ) 


computed from the four-term approximation were 
plotted as a function of m figure 7 for three 


values of eft. The values of cjt were taken as 500, 
312.5, and 240.9, values consistent with practical 
design considerations and within the limits of the 
present idealization for monocoque and multiweb 
wings. As can be seen from figure 7, very little 
effect is noted in the flutter boundary caused by 
variations in cjt. 

Effects of aerodynamic damping. — In both the 
exact analysis and in the Galcrkin approximations 


aerodynamic damping has been neglected. The 
damping term was neglected in the exact analysis 
in order that a simplified solution to the differen- 
tial equation could be found. For the Galerkin 
approximations, however, the inclusion of the 
aerodynamic damping term docs not create any 
difficulties. Thus, the effect of aerodynamic 
damping on the critical flutter parameter has been 
investigated for the limiting conditions of pin- 
ended (k=™) and free-free (Jc= 0) constraints by 
using a Galerkin modal approach. 

A two-term Galerkin approximation (which is 
identical to the three-term approximation at /r=0) 
was employed for the sake of simplicity, since the 
choice of the proper coupling modes has already 
been established. Tt is assumed therefore that 
the critical coupling modes do not change with the 
inclusion of aerodynamic damping and that only 
the magnitude of the critical A is affected. 
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The determinant of the coefficients (which are 
now complex) for the two-term Galerkin approxi- 
mation reduces to 



where B x is now complex and is given by equation 
(13) and w n S and A are the general vibration 

l 

h 

For&=<», the first an Asymmetrical (rigid-body 
torsional, m = 0) and the first symmetrical chord- 
wise (ft — 1) modes should be used in equation (43). 



modes for 




For k 0 the first an tisym metrical mode is zero 
and uncouples from the solution; the second anti- 
symmetrical (m — 1) and the first symmetrical 
chordwise (n= 1) modes should be used in this case. 

If 7? 4 is complex and co is considered real, expand- 
ing and simplifying equation (43) and equating 
the coefficients of both the real and imaginary 
parts to zero yields t lie following expression for A cr : 

A ct ] h « a) 


+ \ V(«£, s ffll, ,i) 2 + 1 2SKi(a 4 at ,t) 2 ] (44) 

where 


2 Wpvil/*- 


(45) 


and 



In order to arrive at numerical results of the 
effects of aerodynamic damping on the flutter 
parameter A in the present investigation, the data 
from the flutter failures of the multiweb-wing 
structures were used. The pertinent parameters 
used in some of these tests are tabulated below: 

5 = 312.5 

(£) Uo187 

—=0.0001453 

P 

7?= 10.6X10" psi 
M=3.0 


The results obtained by using the tabulated data 

indicate that for k— 0 and — 0 (ai i5 = 4.591; 

*=3.927), the effect of aerodynamic damping 
increases \ cr by only 0.3 percent. If the actual 

^ ^ value of 0.0187 is used in place of the 

approximation, =0, the results arc not al- 
tered; the increase in A cr is still only of the order 
of 0.3 percent. 


When fc— oo and — 0 (<*1,5=4.591 ; a QiA = w), 

the effect of aerodynamic damping on A cr is still 
negligible, the value of A cr being increased by only 
0.2 percent. A greater effect would be realized in 
these limiting cases, however, if the cjt ratios 
increased, as can be seen from the density-ratio 
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factor tj in equation (45). Note that the effects 
of damping have been studied only for the limiting 
cases of k^= 0 and k= °° . It is nevertheless believed 
that the conclusions based on these cases are 
applicable for finite values of k except in the 

questionable region of ^ ~^.35, where flutter 

appears to occur’ for negligible air flow. 

APPLICATION OF METHOD 

The previous sections discussed the flutter be- 
havior of an idealized structure. The purpose of 
the present section is to show how these results 
can be used to predict the Rutter boundary of the 
multiweb wings used in the experimental investi- 
gations. In addition, the accuracy of the predicted 
flutter boundary is ascertained by comparison with 
actual flutter data. 

TEST SPECIMENS 

Tn the experimental investigation into the dy- 
namic flutter failures of multi web-wing structures, 
a series of multiweb wings were tested in the pre- 
flight jet at the NASA Wallops Station at a Mach 
number of 2 (see refs. 1 to 3) and in the Langley 
9 _ by 0-foot thermal structures tunnel at a Mach 
number of 3. Specimen identification (see table 
IT) conforms with that of the multiweb-wing 


(MW) series previously discussed in references 1 
to 3. Numbers within parentheses indicate dupli- 
cates of the same specimen. For example, MW- 
2-(2) identifies the second MW-2 model tested. 
These aluminum-alloy (/£= 10. 6X 10° psi) multi- 
web models had 5-percent-thick circular-arc air- 
foil sections with solid leading and trailing edges 
as shown in figure 1 (a). The chord length of 15 
of the IS models tested (models MW-2 and MW -4) 
was 20 inches and the skin thickness f , 0.004 inch; 
thus the ratio cjt for these specimens was 312.5. 
These models had a 20-inch span (approximately) 
with solid root and tip bulkheads but were free of 
internal ribs. (The models designated MW-2 had 
0.2 5-inch tip bulkheads and those designated 
MW-4 had 0.025-inch tip bulkheads.) The pa- 
rameter was fk0bS7 and the wing aspect ratio 
slightly over 1.0. 

The last three multiweb models listed in table IT 
(designated MW-l-(2), MW-23, and MW-24) are 
presented in order to compare the theoretical flut- 
ter boundary with multi web models of different 
chord length and skin thickness from those of 
models MW-2 and MW-4. The ratio cjt of these 
multiweb models, however, is quite similar to that 
of models MW-2 and MW-4 (c/#=320 toc/*=312.5 
for MW-2- (4)). Model MW-l-(2) (from ref. 2) 


TABLE II. — AERODYNAMIC TEST CONDITIONS 


[p =0.000259 slug/eu In.] 


Model 

Mach 

number 

Stagnation 
temporal urn, 

°F 

Velocity, 

r„. fi« 

Air density, 

Pao , slugs/cn ft 

Dynamic 
pressure, q r 
lb/sq in. 

Flutter 

References 

MW-2-(2) 

1 99 

89 

1 , 700 

0. 003010 

30. 40 

Yes 

1 

MW-2- (2) 

L 99 

540 

2, 300 

. 002290 

42. 00 

\ es 

1 

MW-2- (2) 

1. 99 

503 

2, 2G0 

. 002320 

41. 10 

Y es 

1 

MW-2 -(2) 

E 99 

517 

2, 280 

002300 

41. 50 

V es 

1 

MW-2 (2) 

E 99 

474 

2, 230 

. 002380 

40. 90 

Yes 

1 

MW-2- (2) 

1. 99 

to 

Zj\ 

2, 200 

. 002300 

41. 80 

A es 

1 

MW-2- (4) 

1. 99 

551 

2, 314 

. 002220 

41. 44 

V es 

3 

MW-4- (3) 

1. 99 

531 

2, 288 

. 002280 

41. 58 

Y es 


MW-4- (4) 

3. 00 

651 

2, 928 

. 001095 

32. G2 

Yes 

En published data. 

MW 4- (5) 

3. 00 

301 

2, 423 

. 001640 

33. 49 

No 

Fn published data 

MW- 4- (5) 

3. 00 

395 

2, 360 

. 001423 

32. G2 

No 

Unpublished data 

MW-4-(6) 

3. 00 

502 

2, 727 

. 001315 

33. 89 

No 

En published data 

MW-4 -(G) 

3. 00 

G0S 

2, 871 

. 001 183 

33. 89 

No 

I T n published data 

MW-4 -(7) 

3. 00 

C73 

2, 9G1 

. 001140 

34. 71 

Yes 

T T n published data 

MW-4- (9) 

3. 00 

G59 

2, 943 

| • . 000780 

23. 62 

No 

Enpublished data 

MW-l-(2) 

E 99 

574 

2, 350 

. 002120 

40. 60 

Yes 

2 

MW-23 

3. 00 

G47 

2, 923 

. 001140 

33. 82 

No 

En published data 

MW-24 

3. 00 

G5G 

2, 935 

. 001077 

32. 24 

Yes 

En published data 
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has a 40-inrh chord, skin thickness of 0.125, tip 
bulkhead of 1 inch, and a span of 37.50 inches. 
Model MW-23 has the same dimensions as MW- 
l-(2) except for a tip bulkhead of 0.051 inch and 
a span of 43.25 inches. Model MW-24 has a 
00-inch chord, skin thickness of 0.188, tip bulk- 
head of 0.073 inch, and a span of 04.88 inches. 

(ai ilS =4.56S) and the wing aspect ratio is approxi- 
mately 1.0. 

Tabulated in table II are such pertinent param- 
eters of the experimental tests of these models as 
Mach number, stagnation temperature, velocity 
and density of the jet stream, and dynamic pres- 
sure, Whether or not flutter occurred is indicated. 
The published references from which data were 
obtained are listed in case more information is 
desired. 


) 


for the three models is 0.0178 


The ratio 


( 


EXPERIMENTAL RESULTS 

Although all the multiweb wings had 5-percent- 
thick circular-arc cross sections, the curvature of 


these cross sections can be very accurately ap- 
proximated by the parabolic equation h=h 0 ( 1 — f 2 ). 
Furthermore, the ratio c/t of the MW-2 and MW 4 
multiweb wings is 312.5 and thus the theoretical 
flutter boundary presented in figure 5 (For the 
four-term approximation) is the predicted flutter 
boundary for 15 of the experimental models and is 
again shown in figure 8. in figure 8 are also shown 
the experimental flutter points tabulated in table 
TIT for all the MW models. (Note from table III 
and fig. 8 that some of the flutter points fall in 
such close proximity to each other that they arc 
not easily distinguished.) In order to tabulate the 
experimental flutter points the values of the first 
torsional frequency w r and the first symmetrical 
chord wise frequency u LS must be known. The 
first torsional frequency o) T was determined ex- 
perimentally and is listed in the first column of 
table ITT. The constant value of u ltS was deter- 
mined from equation (24) once a l s was found 
from the general symmetrical frequency equation 

(Wen in ref. 4 for any value of ^ to be 4.566. 



Figure 8. Comparison of theoretical flutter boundary with experimental data. c/t = 312.5. 
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TABLE III.— VIBRATION CHARACTERISTICS FOR FLUTTER SPECIMENS 

[p =0.0002. r )9 slng/cii in.] 


Model 


wr, exp. 


in. 5 


k, lb/in. 



1/3 


MW-2-(2) 
MW- 2- (2) 
MW-2-(2) 
M W-2- (2) 
MW-2-(2) 
MW-2-(2) 
MW-2-(4) 
MW- 4— (3) 
MW-4-(4) 
MW- 4- (5) 
M W-4- (5) 
MW-4-(G) 
MW- 4- (6) 
M W-4- (7) 
MW- 4- (9) 
MW-l-(2) 
MW- 23 
MW 24 


139 

139 

139 

139 

139 

139 

147 

149 

145 

142 

142 

143 

143 

144 
144 

74 

76 

54 


98. 135 
98. 135 
98. 135 
98. 135 
98. 135 
98. 135 
98. 135 
98. 135 
98. 135 
98. 135 
98. 135 
98. 135 
98. 135 
98. 135 
98. 135 
1737. 91 
1525. 53 
7798. 21 


92. 76 
92. 70 
92. 76 
92. 76 
92. 76 
92. 76 

103. 77 

111. 39 
105. 49 
101. 17 

101. 17 

102. 60 
102 . 60 

104. 04 
104. 04 
121. 62 

112. 62 
129. 17 


1. 3824 
1. 3824 
1. 3824 
1. 3824 
1. 3824 
1. 3824 
1. 5465 
1. 6600 
1. 5721 
1. 5077 
1. 5077 
1. 5290 
1. 5290 
1. 5505 
1. 5505 

1. 9428 
l 7990 

2. 0468 


0. 2542 
. 2423 
. 2441 
. 2433 
. 2445 
. 2427 
. 2434 
. 2431 
.3111 
. 3084 
. 3111 
. 3072 
. 3072 
. 3047 
. 3465 
. 2393 
. 3002 
. 3058 


The spring constants k were computed from the 
first torsional frequency by the well-known rela- 
tion 

wr*-^ (47) 

Pip 

where K r is the torsional spring constant and I p 
is the polar moment of inertia per unit length in 
the spanwise direction. From equation (40) 


necessitate applying another set ol elastic re- 
straints horizontally, at the leading and trailing 
edges of the wing cross section, fig. 1 (b), since the 
edges are assumed to remain stationary in the 
present analysis.) However, the magnitude of 
the effect of these thermal stresses on the flutter 
behavior of the multiweb wings is not too evident 
from a comparison of the experimental data and 
the theoretically predicted flutter boundary. 


K r -2k(g 

Thus 

2 uSpTj, (4g) 

c 2 

As can be seen from figure S the predicted 
flutter boundary and the experimental flutter 
points agree quite well. Some of the experi- 
mental no-flutter points do fall slightly below the 
flutter boundary in the flutter region, but it is be- 
lieved that this may be due to the omission of 
temperature considerations in the calculated 
flutter boundary. The heated air ol the test fa- 
cilities produces a temperature gradient between 
the skin and webs of the multiweb model inducing 
thermal stresses. No attempt lias been made to 
include these thermal stresses in the present anal- 
ysis. (The inclusion of thermal effects would 


CONCLUDING REMARKS 

The flutter characteristics determined from an 
idealized model have been used successfully to 
predict camber or cross-sectional flutter behavior 
of actual multi web and monocoque wings of low 
aspect ratio. The use of modified piston theory 
for the air loads in the flutter analysis yielded 
flutter results that agreed well with existing ex- 
perimental data. An approximate three- and four- 
term Galerkin procedure utilizing (lie exact modes 
and frequencies of the vibrational analysis of those 
multiweb wings was employed to establish the 
flutter boundary. The correlation between this 
approximate solution and an exact solution for the 
specific case of simply supported end constraints 
indicated that the inclusion of more terms in the 
approximate modal solution was not necessary , 
the three- and four-term approximations being 
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within it few perron! of the exact. This agree- 
ment was also borne out by the rapid convergence 
of the Galerkin procedures for the three- and four- 
term approximations. 

The existence of curvature in a wing cross sec- 
tion was shown to have an appreciable effect on 
the flutter criterion. Once the curvature is pres- 
ent, however, variations in the degree of curvature 
were shown to have little effect on the flutter 
boundary. The effects of aerodynamic damping 
on the flutter boundary were investigated for (he 
limiting cases of £—0 and k= °° and found to be 


negligible. 

The influence of thermal stresses arising in the 
multi web-wing structure due to heating of the 
skins by the hot-air jet have not been considered 
in the present analysis, but a comparison of t lie 
experimental data with the theoretically predicted 
flutter boundary did not indicate the results to 
be unduly affected by thermal stresses. 

Langley Research Center, 

National Aeronautics and Space Administration, 
Langley Station, Hampton, Va., Auguti 23, 1962. 



APPENDIX 


EXACT SOLUTION OF DIFFERENTIAL EQUATION FOR k= AND AERODYNAMIC DAMPING NEGLECTED 


An exact solution to equation (8), as given in 
the text, can be obtained directly by neglecting 
aerodynamic damping. In a manner similar to 
that employed by Hedgepeth in reference 6, an 
exact solution is obtained for the limiting case of 
pin-ended constraints at the leading and 

trailing edges of the cross section. With aero- 
dynamic damping neglected, /? 4 reduces to —a 4 and 
equation (8) can be written as 


10 IV + Aw f — a*w- 


— tiXV 2 
Et 3 


(Al) 


Furthermore, the origin of the coordinate axes is 
moved to the leading edge of the cross section both 
to aid in a later comparison with the infinite- 
aspect-ratio simply supported plate and to simplify 
somewhat the form of the solution to equation 
(Al). Thus 

*=f+l (A2) 

and the expression for the half-depth h can be 
written as 

h=K [l-(*-D*] (A3) 

The associated boundary conditions for the pin- 
ended constraints become 


w(0)=w(2) = 0 


(A4) 


where (for p— 1 to 4) are the roots of the 
auxiliary equation 


T A m — a A = 0 


(AS) 


Since w is considered real, a 4 is always positive and 
A is always positive by definition; hence, m lias 
one possible real positive root and one possible 
real negative root and a pair of complex roots, 
or two pairs of complex roots. 


Let 


7U] — 0T 4> ^ 

m 2 = 0— 0 l 

m*=y+ih f 
m 4 — y — id J 


(A9) 


In order for these roots (eqs. (A9)) to satisfy 
equation (AS), the following equations must be 
true 

m j + m 2 + m 3 + m A =■ 0 (A 1 0) 

m 1 w 2 +m 3 m 4 + (w, + m 2 ) (wi 3 + m 4 ) =0 (Al 1) 

mpnAmy; m 4 ) + m 4 m 3 (m I + m 2 ) — —A (Al2) 

m i m 2 m -pn 4 = — a 4 (A 1 3) 

By substitution of the relations for m Xy m 2} m 3 , and 
jrii given in equations (A9) into equations (A 10) 
to (A 13) the following relationships between 7, <5, <f> } 
and 0 are found: 


a w 

77^ 

u w 

=0 (A5) 

*»2 

y=-0 

(A 14) 

The normal midplane force X is considered con- 

<5 2 = + 2(P 

(A 15) 

stant and equal to the average midplane force as 

20(</> 2 + 5 2 ) = -A 

(A 16) 

given by equation 
simply reduces to 

(15) or, 

since u'(0) ~ w(2) = 0, 

{()- -- 0 2 ) (0 2 + 5 2 ) = -a 4 

(A 17) 


4A 0 FA 


Vvxty 

Jo 


(AO) 


The solution to equation (Al) can be written in 
the form 


By inspection of equations (A14) to (A 17) it can 
be seen that 

0 < 0 \ 

!<#»]> \e\J 


(A 18) 


_ nx/(„c 2 , ^ 
W ~'a'Et* + S 


(A7) 


while from equation (A 17) the possibility of two 
pairs of complex roots is eliminated. Substitu- 
tion of equation (A 15) into equations (A 10) and 
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(A17) yields, after some algebraic manipulation, 

(Al« 

and 

#+\ (A20) 

In equation (A20), 0 2 lias one positive real root 
and two conjugate complex roots. However, 
interest is centered in the case where 6 is real and 
negative (eq. (AIR)); hence Q 2 is real and positive. 
Substitution for <p 2 from equation (A19) into 
equation (A 15) yields 

‘■=**-5 < A 21 ) 

Therefore, for the flutter analysis, equation (A7) 


can be written as 

+e*(C, cosh M+C* sinl) 

or lit 

d -e~ 9 * {C $ cos fy+C A sin 50) (A22) 

Substitution of the general solution (given by eq. 
(A22)) into the pin-ended boundary conditions 
(eqs. (A4) and (AS)) yields the following set of 
equations for the coefficients C x , C and C A : 

0 (.123) 

dfeVCc 29 cosh 24 >-{-C 2 e 2S sinli 2 <f> 

+ C s e~ 2e cos 2 d+C 4 e~ 2e sin 25=0 (A24) 


C x (0 2 + 0 2 ) + G200+ G(0 2 - 5 2 ) - cm = o 


Ci[(0 2 + 0 2 )c 20 cosh 20+2^0c 20 sinh 20]d Ca[2^c 2tf cosh 20 + (0 2J r<f> 2 )e 29 sinli 20] 

+ ft[(0 2 -5 2 )e~ 25 cos 26+268er 2d sin 25] — C r 4 [2^«“ 2 * cos 25- (0 2 — 5 2 )c“ 2 * sin 25] = 0 (A26) 

where X is written as 



(A27) 


For a nontrivial solution the determinant of 
the coefficients must equal zero. Substitution of X 
given by equation (A27) into equations (A23) and 
(A24) and subsequent simplification yields the 
determinant of the coefficients in terms of 6, 0, 5, 

a \ aT1( l A as shown in table I. 


Note that if the half-depth h is equal to zero 
everywhere (h 0 = 0), the parabolic airfoil reduces 
to a flat beam and a in table I becomes zero. This 
modified determinant, when expanded, will yield 
an expression which is identical to that for the 
infinite-aspect-ratio simply supported plate solved 
exactly in reference 6. 
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